Ion temperature profile stiffness: non-linear 
gyrokinetic simulations and comparison with 



experiment 



J. Citrin 1 , C. Bourdelle 2 , J.W. Haverkort 13 , G.M.D. Hogeweij 1 , 
F. Jenko 4 , P. Mantica 5 , M.J. Pueschel 6 , D. Told 4 and 
JET-EFDA contributors* 

JET-EFDA, Culham Science Centre, Abingdon, 0X14 3DB, UK 

1 F0M Institute DIFFER - Dutch Institute for Fundamental Energy Research - 

Association EURATOM-FOM, Nieuwegein, The Netherlands 

2 CEA, IRFM, F-13108 Saint Paul Lez Durance, France 

3 Centrum Wiskunde & Informatica (CWI), PO Box 94079, Amsterdam, The 
Netherlands 

4 Max Planck Institute for Plasma Physics, EURATOM Association, Boltzmannstr. 
2, 85748 Garching, Germany 

5 Istituto di Fisica del Plasma P. Caldirola, Associazione Euratom-ENEA-CNR, 
Milano, Italy 

6 University of Wisconsin-Madison, Madison, Wisconsin 53706, USA 

*See the Appendix of F. Romanelli et al., Proceedings of the 24th IAEA Fusion 
Energy Conference 2012, San Diego, USA 

E-mail: J.Citrin@differ.nl 

Abstract. Recent experimental observations at JET show evidence of reduced ion 
temperature profile stiffness at low magnetic shear (s) in the presence of flow shear. 
Non-linear gyrokinetic simulations are performed, aiming to investigate the physical 
mechanism behind the observations. The sensitivity of profile stiffness to the variations 
of plasma parameters experimentally observed when transitioning to the low-stiffness 
regime is assessed. It is found that non-linear electromagnetic effects, even at low 
f3 e , can significantly reduce the profile stiffness, although not by a degree sufficient to 
explain the experimental observations. The effect of toroidal flow shear itself is not 
predicted by the simulations to lead to a significant reduction in flux due to significant 
parallel gradient velocity destabilisation. For the majority of discharges studied, the 
simulated and experimental ion heat flux values do agree within reasonable variations of 
input parameters around the experimental uncertainties. However, no such reasonable 
agreement is obtained for the discharge with the highest logarithmic ion temperature 
gradient. The simulated stiffness level is thus higher than observed for this regime, 
when assuming pure toroidal rotation. 
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1. Introduction 

It is well established that one of the primary limitations of tokamak energy 
confinement is ion-Larmor-radius-scale turbulent transport driven by background 
pressure gradients pQ. The ion-temperature-gradient (ITG) instability in particular has 
been long identified as a ubiquitous unstable mode in tokamak plasmas (21 [3j HJ [5] , and is 
primarily responsible for ion heat losses. The instability saturates in a non-linear state 
in conjunction with non-linearly excited zonal-flows, forming a self-organised turbulent 
system which sets the transport fluxes [6]. 

In addition to self-organised zonal-flows, the application of external flow shear 
is predicted to suppress turbulence through two broad mechanisms: decorrelation of 
the turbulent structures in the non-linear phase once the shearing rate is comparable 
with or exceeds the inverse non-linear autocorrelation time; suppression of the driving 
linear modes by continuously shifting the mode from the most unstable spatial 
scale to nearby, more stable spatial scales [7J E\- Flow shear has been observed 
experimentally to be correlated with tokamak transport barriers [H [TO], [11] . Non-linear 
gyrofluid simulations with adiabatic electrons have predicted turbulence quenching 
above ^E/lmax = 1±0.3 [El US], where for purely toroidal rotation the normalised ExB 
shear rate j E ~^/(^), and ^ ma x is the maximum linear growth rate in the absence of 
rotation. Later gyrokinetic simulations, including cases with kinetic electrons, predicted 
that quenching occurs at somewhat higher (but similar) flow shear compared to the 
gyrofluid simulations, at "fE/jmax — 2 ±0.5 [131 HSJ LTE] - According to results in Ref . [T5] . 
this quench behaviour is independent of the adiabatic or kinetic electron assumption. 
However, when including kinetic electrons, the trapped electron drive tends to raise the 
instability growth rates even for ITG turbulence. Therefore, when including kinetic 
electrons (which is more realistic) the resultant higher ^ max necessitates a higher value 
of 7s compared with the adiabatic electron case to reach a similar •jE/lmax ratio and 
quench the turbulence. The seeming robustness of the transport quench has motivated 
formulations of effective growth rate reduction due to the flow shear in the mixing length 
rule of quasilinear transport models such as GLF23 and TGLF [T7J [18] . 

The abovementioned quench results were obtained in simulations which did not 
include a self-consistent parallel velocity gradient (PVG) in the system, which can be 
destabilising [191 120] ■ When PVG destabilisation is included, simulations have shown 
that they can limit the transport quench [TSJ EH [15] . It is thus important to incorporate 
the effects of PVG destabilisation in transport models when comparing modelling 
predictions with experiments that have significant flow shear. For pure toroidal rotation, 
the degree of the PVG destabilisation depends on the magnetic geometry since then 
7 P = ^)Ei where 7 P is the PVG shear rate, and e=r/R. 

Experimental results of flow shear stabilisation of core turbulence are not 
systematic. While improved core confinement in hybrid scenarios due to flow shear has 
been observed at DIII-D [22], a counterexample is evident from experiments involving 
JET standard H-modes which showed a lack of variation in core-confinement even when 
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significantly varying the rotation profiles and maintaining constant total power [23J. 

It has been recently observed in dedicated experiments at JET that ion temperature 
profile stiffness can be reduced at low normalised radii (r/a < 0.5), in disagreement 
with non-linear ITG turbulence modelling [24"| |25| 126] . This has been hypothesised 
to be related to the correlation between low magnetic shear (s) and increased flow 
shear in the low stiffness discharges. The term 'stiffness' is defined here as the local 
gradient of the gyro-Bohm normalised ion heat flux with respect to R/Lti (normalised 
inverse gradient length). The observations concentrated on p = 0.33 and p = 0.64 
(where p is the normalised toroidal flux coordinate). At p = 0.33, the stiffness it 
observed to transit from high to low when the flow shear was increased. However, at 
p = 0.64, stiffness is observed to be high irrespective of flow shear. A previous non- 
linear gyrokinetic study based on the recent JET discharges at p = 0.33, as detailed in 
Ref . [25] . reported only a ITG threshold shift with rotation, as opposed to a decrease in 
stiffness as observed. Additional observations made in Refs. [2U [251 126] pertinent to this 
work are as follows: at low rotation at p = 0.33, the observed stiffness level is higher 
than the gyrokinetic simulation predictions; furthermore, the observed ITG threshold 
is lower than the non-linear gyrokinetic prediction, questioning the manifestation of the 
Dimits shift [27] predicted by non-linear simulations. 

In this paper, we extend this previous work and aim to investigate whether 
the experimental observations can be understood through gyrokinetic modelling. 
Understanding these effects could allow the identification of a potential actuator for 
core T{ control. For the analysis, linear and non-linear simulations are carried out with 
the Gene code [28J . Four JET discharges (with the previous carbon wall) were selected: 
70084, 66130, 66404, and 73221. Discharges 66130 and 66404 are situated on the 
'high-rotation, low-stiffness branch' at p = 0.33 seen in Fig.l in Ref. [25], and partially 
reproduced here for convenience in FigJTJ Discharge 70084 is a low flux, low rotation 
discharge selected to provide a data point near the turbulence threshold. Discharge 
73221 is a high flux, low rotation discharge situated on the 'low-rotation, high- stiffness 
branch' at p — 0.33, as shown in FigHJ The specific questions which we investigate are 
the following: 

(1) Is the experimentally observed stiffness reduction for the high-rotation discharges 
at p = 0.33 consistent with gyrokinetic non-linear simulation predictions? Which 
plasma parameters have the highest impact on the stiffness level for ITG turbulence? 
Is there sufficient leeway in the plasma parameters due to uncertainties such that the 
experimental observations and non-linear simulation predictions can be reconciled? 

(2) Is the lack of experimentally observed stiffness reduction for the high-rotation 
discharges at p = 0.64 consistent with gyrokinetic non-linear simulation predictions? 

(3) Can the experimentally extrapolated turbulence threshold be reconciled with 
the non-linear turbulence threshold including the Dimits shift, given reasonable 
uncertainties in the plasma parameters? 

(4) Can the seeming high stiffness in the 'low-rotation, high-stiffness' branch at p = 0.33 
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be reconciled with the non-linear simulations, given reasonable uncertainties in 
plasma parameters? 

The discharges were reanalysed with the CRONOS suite of integrated modelling 
codes [29] to identify any differences in parameters apart from rotation and R/Lt% 
within the chosen discharge set - such as T e /T,, R/L n , (3 e , and fast particle content - 
which may lead to the observed differences in ion heat flux and R/Lxi- The sensitivity 
of the ion heat flux and stiffness to each of these parameters was tested with Gene in 
dedicated R/L^i scans. Finally, complete simulations - i.e. collisional, electromagnetic, 
multi-species, and with rotation - were carried out at both p = 0.33 and p = 0.64. 

The rest of the paper is organised as follows. In section [2] the Gene gyrokinetic 
code is briefly reviewed, as are the base parameters of the simulated discharges. In 
section [3] the results at p — 0.33 are shown. In section @] the results at p — 0.64 are 
shown. Conclusions are presented in section |5j 

p=0.33 p=0.64 




02468 10 02468 10 



Figure 1. Partial reproduction of data presented in Ref. |25j displaying the separation 
between high and low stiffness regimes at p = 0.33 (a) for discharges with low 
and high rotation respectively. At p = 0.64 (b) no significant separation of the 
stiffness behaviour is evident. The heat fluxes are in gyroBohm normalised units, 
qGB — T?- 5 n i m ( -- 5 /e 2 B 2 R 2 . The specific discharges studied in this paper have been 
circled. 



2. GENE simulations and discharge parameters 

Gene solves the gyrokinetic Vlasov equation, coupled self-consistently to Maxwell's 
equations, within a 5f formulation [30]. Computational efficiency is gained by solving 
in field line coordinates, x is the radial coordinate, z is the (poloidal) coordinate along 
the field line, and y is the binormal coordinate. Both an analytical circular geometry 
model (derived in Ref.|31j) as well as numerical geometry were used in this work. The 
circular geometry model avoids the order e = a/R inconsistency present in the often 
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Table 1. Discharge dimensional parameters. The values are averaged between 9.5-10.5 
s for discharges 70084 and 66130, between 6.5-7.5 s for discharge 66404, and between 
7-8 s for discharge 73221. Quoted errors are statistical, and do not include possible 
systematic errors. 



Shot no. ©location 


B [T] 


I P [MA] 


Ti [keV] 


T e [keV] 


n e [10 19 m" 3 ] 


70084@p = 0.33 


3.5 


1.8 


2.01 ±0.02 


2.16 ±0.1 


2.6 ±0.2 


66130@p = 0.33 


3.1 


1.5 


2.58 ±0.04 


3.3 ±0.3 


3.37 ±0.24 


66404@p = 0.33 


3.5 


1.8 


3.1 ±0.13 


3.61 ± 0.08 


2.3 ±0.1 


73221@p = 0.33 


3.5 


1.8 


1.84 ±0.04 


2.44 ± 0.03 


2.45 ±0.16 


70084% = 0.64 


3.5 


1.8 


1.08 ±0.02 


1.27 ±0.05 


2.3 ±0.3 


66130@p = 0.64 


3.1 


1.5 


1.38 ±0.03 


1.5 ±0.24 


2.8 ±0.3 


66404% = 0.64 


3.5 


1.8 


1.34 ± 0.08 


1.66 ±0.14 


1.51 ±0.07 



applied s — a model, but does not include a Shafranov shift. For the numerical geometry, 
the finesse code was used to solve the extended Grad- Shafranov equation including 
toroidal rotation [32] ■ All simulations carried out were local, justified since l/p*~500 
for the range of plasma parameters studied here [331 [M]. Both linear and non-linear 
simulations were performed. In the linear mode, an eigenvalue solver was used to 
compute multiple modes for each point in parameter space [35| 136] . In the presence 
of rotation, when no time-independent eigenmodes can form, a complementary initial 
value solver was used. 

Four discharges from the data-set presented in Ref.[25] were analyzed at p = 0.33 
and p = 0.64, where p is the normalised toroidal flux coordinate. The discharges are 
70084, 66130, 66404, and 73221. Discharge 70084 corresponds to a representative low 
rotation, low flux discharge. 66130 and 66404 are discharges further up on the 'high 
rotation, decreased stiffness' curve as seen in FigJTJ 73221 is a high flux, low rotation 
discharge situated on the 'low- rotation, high-stiffness branch' at p = 0.33, as shown in 
Fig{TJ The kinetic profiles of the four discharges were spline fitted and interpretative 
runs were carried out with the CRONOS integrated modelling suite of codes [29J for 
the equilibrium calculations and g-profile calculations. The kinetic profiles were then 
averaged over 1 s centered around 10/10/7/7.5 s respectively for calculations of the 
gradient lengths and other quantities such as /3 e . The parameters are shown tables 
[TK2J Discharge 73221 was only analysed at p = 0.33, for the investigation of the 
seemingly high stiffness of the low rotation branch. The (Z e ff) values correspond to 
Bremsstrahlung measurements. Since the precise Z e ff profiles are not known, the 
sensitivity of the transport predictions to the range of reasonable Z e ff at p = 0.33 
is explored in section [3771 is* is the normalised collisionality: u*=v e i , with e = a/R 

and v te = \f^f- Note that the data presented in table |2] was processed separately and 
independently from the values quoted in Ref.[2U [25] and shown in FigfTJ The R/Lt{ 
values in table |2] and FigJT] agree within error bars. 

The agreement between the g-profiles obtained by CRONOS interpretative 
simulations and the measured g-profiles is satisfactory, as seen in Figj2j The average 
discrepancy between the interpretative and measured g-profile values at p = 0.33 
and 0.64 is ~ 10%, within the estimated uncertainty of the g-profile measurements. 
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Table 2. Discharge dimensionless parameters. The s and q values are calculated by 
CRONOS interpretative simulations, assuming neoclassical diffusion. The values are 
averaged between 9.5-10.5 s for discharges 70084 and 66130, between 6.5-7.5 s for 
discharge 66404, and between 7-8 s for discharge 73221. 



Shot no. ©location 


8 


3 


Te/Ti 


R/Irn 


R/LTe 




ft [%] 




(Zeff) 


Mlvtar/cs] 


70084@p = 0.33 


0.7 


1.7 


1.08 ±0.04 


3.5 ±0.5 


3.8 ±0.6 


1.4 ±0.4 


0.19 ±0.01 


0.07 


2.2 ±0.1 


0.09 


66130@p = 0.33 


0.7 


1.8 


1.25 ±0.13 


6 ±0.4 


6.5 ± 1 


2.4 ±1 


0.46 ±0.09 


0.04 


1.8 ±0.1 


0.31 


66404®/!) = 0.33 


0.4 


1.8 


1.14 ±0.06 


8.6 ±0.9 


5.5 ±0.8 


3.8 ±0.4 


0.35 ±0.07 


0.02 


2.2 ± 0.1 


0.19 


73221@p = 0.33 


0.7 


1.5 


1.33 ±0.02 


3.8 ±0.4 


5.4 ±0.2 


2.8 ±0.3 


0.2 ±0.02 


0.055 


2.2 ±0.1 


0.07 


70084Cojp = 0.64 


1.3 


3 


1.18 ±0.05 


7.2 ± 0.2 


6.4 ± 1 


1.8 ±0.8 


0.096 ±0.01 


0.16 


2.2 ± 0.1 


0.03 


66130@p = 0.64 


1.5 


3.5 


1.1 ±0.2 


6.8 ±0.3 


8.5 ±3 


1.8 ± 1.4 


0.18 ±0.04 


0.1 


1.8 ±0.1 


0.23 


66404@p = 0.64 


1.4 


2.9 


1.23 ±0.13 


6.9 ±0.4 


10 ± 1.6 


2.1 ±0.9 


0.08 ±0.01 


0.05 


2.2 ±0.1 


0.15 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



P P P P 

Figure 2. Comparison between CRONOS interpretative simulation q-profiles and 
experimental q-profiles. The profiles are averaged between 9.5-10.5 s for discharges 
70084 and 66130, between 6.5-7.5 s for discharge 66404, and between 7-8 s for discharge 
73221. 



The experimental g-profiles were obtained by EFIT constrained by either Faraday 
rotation measurements (discharges 70084 and 73221) or motional Stark effect (MSE) 
measurements (discharges 66130 and 66404). 

In the Gene simulations, typical grid parameters were as follows: perpendicular box 
sizes [L x ,L y ] = [170, 125] in units of p s =c s /Q ci = ^T e /mi/ (eB/nii), perpendicular grid 
discretisations [n x ,n y ] = [192,48], 24 point discretisation in the parallel direction, 32 
points in the parallel velocity direction, and 8 magnetic moments. Extensive convergence 
tests were carried out for representative simulations throughout the parameter space 
spanned in this work. The lack of convergence of the heat fluxes with increasing n y 
as reported for gyro [37J simulations of discharge 70084 in Ref. [25] - associated with 
increasing zonal flows - was not encountered here. In our cases the convergence with 
n y was well behaved. The difference may stem from the different treatment of the 
radial boundary conditions in the Gene and GYRO simulations. Further investigation 
is necessary to ascertain this. The heat fluxes shown in the following sections are in 
gyroBohm normalised units, qcB — T 2 - 5 nimf 5 /e 2 B 2 R 2 . k y is in units of l/p s . These 
heat fluxes correspond to time averaged values over the saturated state of the Gene 
simulations. The statistical flux variations due to intermittency are for clarity not 
explicitly shown as error bars. This variation is typically 5 — 10% for our parameters. 
7 and are in units of c s /R where c s = JT e /mi. All rotation is considered to be 
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purely toroidal unless specifically mentioned otherwise. For the low and high rotation 
discharges it was assumed that j E = 0.1 and 0.3 respectively, at both p = 0.33 and 
p = 0.64. These are representative 'Je values for the low and high stiffness discharges 
from the dataset in Ref . [25] . 

3. Simulations at p = 0.33 

In this section, we isolate the effect of various parameters on ion profile stiffness 
and critical threshold, at p = 0.33 (where the transition to low stiffness at high 
rotation was observed). These parameters are: q, s, rotation, effect of rotation on the 
magnetohydrodynamic (MHD) equilibrium, fast particle content, Rj L n , (3 e , and Z e ff. 
We then proceed to realistic simulations of 70084, 66130, 66404 and 73221 for a full 
comparison between the gyrokinetic predictions and experimental ion heat fluxes. These 
simulations simultaneously include: numerical geometry, collisions, electromagnetic 
effects, Z e jf (by including a carbon species), and realistic T e /Tj. 

3.1. Stiffness and threshold sensitivity to q and s 

While the linear ITG turbulence threshold increases with s/q [38], the stiffness (i.e. 
the rate of change of the gyro-Bohm normalised ion heat flux with respect to R/Lti) 
decreases in non-linear ITG simulations with both decreasing s and decreasing q. For 
decreasing s, the reduced stiffness has been shown to be correlated with increased 
coupling to zonal flows [39]. For decreasing q, this is due to an decreased downshift 
(compared with the peak in the linear spectrum) in the peak wavenumber of the 
turbulence spectrum, indicating decreased correlation lengths [UH HU 02]. These 
sensitivities are shown in FigJSJ For example, we can see that for both the s/q = 0.6/1.3 
and s/q = 1/2 cases the turbulent threshold is similar while the stiffness is lower for the 
s/q = 0.6/1.3 case. 

We will deliberately make an optimistic assumption that s/q = 0.2/1.3 throughout 
all the subsequent parameter scans carried out at p = 0.33 in this work. For the 
numerical geometry cases, this was done by modifying the current profile input into 
CRONOS such that at p = 0.33 values of s/q — 0.2/1.3 were obtained following the 
solution of the Grad-Shafranov equation. The choice of assuming s/q = 0.2/1.3 is to 
ensure that we are in a 'low-s regime', which has been hypothesised to be an important 
factor in the stiffness reduction, based on the observed correlation between low stiffness 
and low-s throughout the data set in Ref. [25] . An added advantage is that this choice is 
optimistic regarding stiffness reduction. Thus, if the non-linear gyrokinetic simulations 
do not display the same degree of low stiffness as experimentally observed even with 
s/q = 0.2/1.3 - as will indeed be shown - then this choice ensures that uncertainties in 
the g-profile cannot be invoked to suggest that the s and q values used were not low 
enough to obtain reduced stiffness. For the Gene simulations at p = 0.64 shown in 
section HI the CRONOS calculated q and s values were taken for each discharge. 
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Figure 3. Non-linear electrostatic collisionless Gene R/Lti scans for various levels 
of s and g-profile with circular geometry at p — 0.33. R/Lxe — 5, and R/L n = 1.1. 



The discussion of the sensitivity of the linear threshold to q brings us to an 
important point. In Refs.JHJ [25], it was pointed out that the measured turbulence 
threshold of the low-rotation discharges in the data set were lower than the predicted 
non-linearly upshifted (Dimits shift) [27J thresholds. These thresholds were predicted 
by non-linear simulations based on discharge 70084 performed with the GS2 non-linear 
gyrokinetic code [IB]. This result questioned the Dimits shift paradigm. The q value 
used for these previous simulations was q = 1.3, based on the processed data at the 
time. However, the data processing methodology for obtaining g-profiles using Faraday 
rotation constraints at JET [H] has since been improved, leading to a revision of the 
measured g-profile value to q = 1.7 at p = 0.33 for t ~ 10 s for discharge 70084. The 
impact of this difference in q on the linear and non-linear thresholds as predicted by the 
gyrokinetic codes is significant. This is shown in FigJH In the figure, the GS2 predicted 
ion heat fluxes for the s/q = 0.6/1.3 case (as shown in Ref. |24j) is compared with the 
analogous Gene simulations. The agreement between the codes is good, apart from 
the zone near the threshold. This difference is likely to be due to the different methods 
used to calculate the geometry: analytical circular in Gene, and s — a geometry in GS2. 
However, the non-linear threshold for s/q = 0.6/1.3 in both codes is approximately 
R/Lti ~ 4.5, above the experimental threshold from Ref. [25]. These curves can then 
be compared with the R/Lxi scan (carried out with Gene) with the revised values 
s/q = 0.7/1.7. In this case, the linear threshold is R/Lti = 2.7, and the non-linear 
threshold following the Dimits shift is at R/Lti ~ 3.5—4, in much better agreement with 
the experimental data. Consistency of the s/q = 0.7/1.7 values with both the revised 
experimental g-profile and CRONOS simulations is thus suggestive that the Dimits shift 
paradigm is in fact now supported by the experimental observations. However, the high 
sensitivity of the turbulence thresholds to the precise s and q values leads us to a more 
conservative conclusion that no firm statement is justified regarding the consistency of 



Ion temp, profile stiffness: non-lin. gyrokinetic simu. and comp. with exp. 



9 



the experimental data with the non-linear Dimits shift. The various values of s and q 
used in the R/L Ti scans in FigJH (including the s/q = 0.2/1.3 values subsequently used 
in this section) can be seen to constitute a sensitivity test of the 'reasonable' range of q 
and s in lieu of rigorous error bars. The one clear conclusion from this sensitivity scan, 
is that there is no clear disagreement between the experimental data and the non-linear 
threshold upshifted due to the Dimits shift. 




Figure 4. Comparison between non-linear electrostatic collisionless Gene and GS2 
R/Lxi scans with the low rotation data from the data-set in Ref.[3S] for various levels 
of s and q. The Gene runs are with circular geometry at p — 0.33, the GS2 runs with 
s — a geometry. R/Ltc — 5, and R/L n = 1.1. 

While the non-linear turbulence threshold extrapolated from the s/q = 0.7/1.7 
curve in FigJD matches the experimental threshold, the simulated stiffness level is 
seemingly lower than the experimental trend. The possibility that this discrepancy 
can be explained by the differences in T e /Tj between the low flux and high flux points in 
the low rotation branch - which impact the critical threshold - is explored in the more 
comprehensive simulations shown in section 13.7.11 

3.2. Stiffness sensitivity to rotation 

In this subsection we isolate the effect of rotation on stiffness, assuming pure toroidal 
rotation. This assumption is justified for JET discharges with significant NBI. 
Collisionless, electrostatic simulations based on 70084 parameters (assuming s/q = 
0.2/1.3) are carried out, applying analytical circular geometry [31]. The predicted 
gyroBohm normalised ion heat fluxes from the R/Lti scans are shown in Fig El The 
sensitivity to 7^ is examined. Even for 7^ = 0.6, double the highest level of flow shear 
achieved in the reference data set from Ref. [25], the simulated level of reduced stiffness 
is significantly less than the experimental observation, as seen by the direct comparison 
with the reference data. However, interesting effects related to the competition between 
stabilising ExB shear and destabilising parallel velocity gradient (PVG) modes - 
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particularly in the vicinity of the threshold - are observed. At low R/Lt%, the PVG 
destabilisation can dominate over the ITG turbulence, reducing stiffness in that region 
of parameter space. Due to the PVG destabilisation, the fluxes do not continue to 
decrease towards the ITG instability thresholds. This is seen in Fig|5k by examining the 
various curves at fixed R/Lt%- At low R/Ltu the fluxes rise with 7^ due to PVG drive. 
However at higher R/L^i, the fluxes decrease with R/Lti due to the ITG stabilisation 
by perpendicular ExB flow shear dominating over the PVG destabilisation. In Figj5b 
the parallel velocity gradients were artificially removed from the system, and the picture 
reverts to a threshold shift. Note that particularly for the (red) je = 0.3 and (black) 
je = 0.6 curves the apparent reduced slope near threshold is not necessarily indicative 
of reduced stiffness in that regime, since the actual effective non-linear threshold may 
lie between the precise values of the R/Lt% values chosen for the simulations. 

For pure toroidal rotation, the relative importance of PVG destabilisation versus 
ExB stabilisation is sensitive to the geometric parameter q/e (where e=r/R) |45j. As 
q/e increases, the field lines are increasingly projected onto the toroidal direction. In 
FiglHl a q/e scan is carried out by varying e in the various R/Lt% scans. Simulations with 
e = 0.11,0.15 assuming circular geometry were performed, as well as an (e)=(r)/R = 
0.13 case from the flux surface averaged minor radius at p = 0.33 using numerical 
geometry from the HELENA [36] equilibrium in the CRONOS simulation of discharge 
70084. The R/Lt% values in the plots corresponding to numerical geometry are defined 
here with respect to the averaged midplane minor radius. The relative strength of the 
PVG destabilisation is seen to weaken as expected with decreasing q/e, until an almost 
pure threshold shift case is reached with q/e = 8.7. 

The interplay between PVG destabilisation and ExB stabilisation demands that 
PVG modes are correctly accounted for in reduced transport models - such as in 
gyrokinetic or gyrofluid based quasilinear models. It is insufficient in such formulations 
to include a growth rate quench model due to Ex B stabilisation without simultaneous 
self-consistent and validated modelling of the PVG modes. Correct modelling near 
the turbulent thresholds is particularly critical for high temperature tokamaks, such as 
ITER. This is because the normalised fluxes are expected to be in the vicinity of the 
turbulence thresholds due to the T^ 2 normalisation dependence. 

Finally, we note that the observed Dimits shift [27] in these cases is only 
A(R/L Ti ) rs 0.5, or alternatively ~ 15%. These values are similar to the 

A(R/L Ti ) ^ shifts observed in previous realistic simulations 07] . The linear 

-ft/ ^Tcrit 1 s 

threshold was calculated by extrapolation to zero-growth-rate of linear R/Lt% scans 
with Gene. The linear threshold in the numerical geometry case is nearly identical to 
the circular geometry case. 

In summary, the Gene simulations do not predict a significant reduction in stiffness 
due to flow shear, even with our deliberate choice of s = 0.2. However, the reduction 
in flux due to the flow shear is not negligible. As suggested by FigEb and as shown 
in section 13.71 a significant reduction of flux due to flow shear is only seen when both 
the effect of PVG destabilisation is artificially reduced, and 7^ is increased beyond the 
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Figure 5. Non-linear Gene R/Lti scans based on 70084 parameters at p = 0.33 
(q/e = 11.8 for circular geometry) and various levels of je [c s /R]- Runs including 
PVG destabilisation are shown in (a). Runs ignoring PVG destabilisation are seen in 
(b). All runs were electrostatic, collisionless, and with circular geometry. The results 
are compared with the low stiffness data at p = 0.33 from Ref.[25 . 



Circular geometry, q/e=1 1 .8 Numerical geometry, q/<e>=1 Circular geometry, q/e=8.7 




23456789 10 23456789 10 23456789 10 

R/L R/L R/L 

Tl Tl Tl 

Figure 6. q/e sensitivity of the PVG destabilisation as seen in R/Lti scans of ion 
heat flux. As q/e is progressively raised, the je induced stabilisation can not only be 
reduced but can even be reversed in the region of the instability threshold. Runs were 
electrostatic, collisionless, and with circular geometry. 

experimental values expected from the toroidal flow shear. 
3.3. Effect of rotation on the equilibrium 

The effect of rotation on stiffness through the impact of the centrifugal force on the 
plasma equilibrium was examined. An extended Grad-Shafranov equation including 
toroidal rotation was solved with the finesse code [32] . using the 70084 pressure and F 
profiles as input, where F=B tor R. For the rotation profiles, scaled variants from 66404 
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were used such that static (7^ — 0), Je — 0.3 and 7^ = 0.6 cases were studied. All 
values correspond to p = 0.33. The different equilibria are seen in FigJTl The sensitivity 
of the equilibria to these levels of rotation are found to be small, as expected due to 
the Mach number squared scaling of the 'rotation pressure'. Only a 10% increase in 
the Shafranov shift was observed between the static and 7^ = 0.6 case. The non-linear 
predicted flux sensitivity to this different Shafranov shift is also minimal, with only a 
6% decrease in ion heat flux when the 7^ = 0.3 equilibrium is used compared with the 
static equilibrium for a run with R/Lt{ = 6.9. We can thus conclude that the effect 
of rotation on the equilibrium itself can only play a minor role in setting the profile 
stiffness. 
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Figure 7. Flux surfaces in the vicinity of the magnetic axis from a solution of the 
generalised Grad-Shafranov equation using the kinetic profiles of 70084 and scaled 
rotation profiles from 66404. Three cases are shown: static (red solid curves), 7^ = 0.3 
(blue dashed curves) and je — 0.6 (black dashed-dotted curves). 



3.4- Inclusion of fast particles 

The discharges studied have relatively low density, and thus for the higher rotation cases 
where significant NBI is employed, it is possible that a large fast ion fraction is present. 
This is predicted to reduce the turbulent drive, since the fast ion species behaves mostly 
adiabatically due to the higher temperatures, and the density of the main ion species is 
reduced through dilution. In ASDEX Upgrade strong evidence has pointed to a fast ion 
dilution mechanism for ITB formation at low density jH]. Monte Carlo simulations of 
the NBI injection and subsequent fast ion slowing down were carried out for discharge 
66404 with nemo/spot jl9] within the CRONOS integrated modelling framework. An 
average fast particle energy (f» 35keV) at p = 0.33 was calculated. A linear Gene scan 
of fast particle densities (relative to n e ) can be seen in FigJHJ In the Gene simulations, 
the fast particle temperature was approximated to the average fast particle energy value. 
The scan is carried out for various k y values in FigJHk, assuming R/Lxfast = 0. The 
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R/L T f ast sensitivity is examined in FigJHb at k y = 0.4. A suppression of the growth 
rates is observed with increasing nf ast /n e . However, for our parameters the fast ion 
fraction is predicted by nemo/spot to be only ~ 10%, and thus insufficient to strongly 
affect the growth rates. While a fast ion fraction has been proposed to be responsible 
for some mismatch between gyrokinetic simulations and experiments [50], in our case it 
seems that the fast ion fraction may be too low to play a significant role. 




Figure 8. Linear fast particle density scans on 66404 parameters at p — 0.33 at 
various values of k y . Runs were electromagnetic, with collisions, and with numerical 
geometry. 



While the nf ast /n e ratio may be low, the fast particle pressure due to NBI and 
ICRH can still be a significant fraction of the total pressure. For the NBI driven 
fast ions in discharge 66404 at p = 0.33, Pf ast /Ptot = 0.21 according to nemo/spot 
modelling (averaging between 6.5-7.5 s). This leads to an increased Shafranov shift 
which, particularly at low magnetic shear, can lead to turbulence stabilisation [51J. 
The effect of the ICRH fast ions on the Shafranov shift was seen to be negligible for 
this discharge. The increased Shafranov shift can be seen in Figj9l where the flux 
surfaces for the low power discharge 70084, and the high power discharge 66404 (with 
and without the inclusion of the NBI fast particle pressure) are compared. The fast 
particle contribution to the 66404 Shafranov shift is significant. For 70084, the Shafranov 
shift is ~ 7.5 cm. For 66404 with the thermal pressure contribution only, the Shafranov 
shift is ~ 8.8 cm. For 66404 with the total pressure (including fast particles), the total 
Shafranov shift is ~ 13 cm. The impact of this difference on the predicted fluxes was 
investigated through dedicated non-linear simulations. The impact was observed to be 
not negligible but also not a dominating factor. A flux reduction of 15% was observed 
in the non-linear simulations with R/L-n = 8 when substituting the numerical geometry 
from 70084 with that of 66404 (i.e. with the fast particle content), as seen in FigJTOl 
We can thus conclude that the effect of fast particles, through both dilution and an 
increased Shafranov shift, can reduce the ion heat flux by ~ 25% when transitioning to 
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the parameters characterising the 'low-stiffness' discharges. This value is not negligible, 
but cannot be the sole explanation for the reduction in stiffness observed. 




R[m] 



Figure 9. Flux surfaces in the vicinity of the magnetic axis for discharge 70084 (red 
solid curves), 66404 without fast particle pressure, (blue dashed curves) and 66404 
with the inclusion of fast particle pressure (black dashed-dotted curves), x = y = 
corresponds to the geometric axis. 




Figure 10. Flux reduction as a function of Shafranov shift normalised to the major 
radius for the three equilibria presented in Fig|9l 



3.5. Impact of R/L n on the stiffness level 

In the limited experimental data set studied, there is a wide variation in R/L n , from 1.4 
in the 70084 case to 3.8 in the 66404 case (which corresponds to the highest R/Lti in the 
data set). The sensitivity of the turbulence to the R/L n value was thus examined. In 
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particular, the possibility that non-linear ITG-TEM (trapped electron mode) interplay 
takes place which can reduce the level of turbulence and thus the stiffness, as reported 
in Ref. [52], was investigated. In FigJTTl these linear scans are shown. For R/L n = 1, the 
dominant mode propagates in the ion diamagnetic direction (ITG mode). However, for 
R/L n = 3.8 the mode at low R/Lxi propagates in the electron diamagnetic direction. 
This is most probably a density gradient driven TEM mode, which is stabilised by R/L^i 
(which would correspond to low stiffness) until it switches to an ITG mode at Rj L Ti w 5. 
At that point we would expect turbulence stabilisation according to Ref. [52] • However, 
for higher R/Lxi the growth-rate stiffness is similar to the R/L n = 1 case, as a pure 
ITG regime is reached. For R/L n = 5 the TEM-dominated regime is maintained for 
much higher R/L^i- However, the highest experimental R/L n in the data set of Ref. [25] 
is R/L n 4. Furthermore, at the experimental high R/Lt% values the transport is ITG 
dominated and stiff even for R/L n = 5. Thus it is unlikely that R/L n is responsible 
for reduced profile stiffness. Furthermore, even if the stiffness is low, the actual growth 
rates themselves are high, and we may expect a high degree of transport. 
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Figure 11. Linear R/Lti scans based on 66404 parameters at p — 0.33 with varying 
R/L n . Growth rates are shown in (a), and frequencies in (b). Runs were electrostatic, 
with collisions, and with circular geometry. 

These results are maintained in the non-linear scans, seen in FigJT2l While at lower 
R/Lxi stiffness is indeed reduced in the TEM regime for the high R/L n case, at higher 
R/Lxi values the difference in stiffness between the R/L n = 1 and R/L n = 3.8 cases 
becomes negligible. We can conclude that the variance of R/L n in the data set is not 
responsible for the observed difference in stiffness. 



3. 6. Impact of (3 e on the stiffness level 

In this subsection the sensitivity of the stiffness on electromagnetic effects - which arise 
for f3 e > - is examined. The simulations carried out take discharge 66404 parameters 
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Figure 12. Non-linear R/Lri scans based on 66404 parameters at p — 0.33 with 
varying R/L n . Runs were electrostatic, with collisions, and with circular geometry. 



as a reference. Linear (at k y = 0.4) and non- linear /3 e scans are shown in FigflBl From 
the linear scans, it is clear that the range of experimental (3 e values (0 — 0.5%) are 
significantly below the kinetic ballooning mode (KBM) thresholds, characterised in the 
plot by the sharp upturn in growth rates at (3 e ~ 1.5 — 2.4%; this finding is expected 
to carry over to the non-linear physics [53] . Below the KBM threshold, j3 e stabilises the 
ITG mode [54J. For our parameters, this leads to a growth rate reduction of ~ 25% 
at {3 e = 0.5%. This is at the upper range of our experimental (3 e values. The 25% 
growth rate stabilisation factor is not exceeded when repeating the linear simulations 
for k y = 0.1 — 0.3. Interestingly, it appears that the linear ITG mode is stabilised at 
lower and lower (3 e as R/L Ti is increased. For practical purposes, however, this should 
be of little benefit as the strong background gradient dependence of the (3kbm threshold 
makes any ITG-stable regions irrelevant, since the KBM threshold significantly decreases 
with R/LTi- 

A striking observation is that the non-linear (3 e ITG stabilisation significantly 
exceeds the linear stabilisation. This is consistent with Gene results reported in 
Refs. [531 EH], as well as, to some degree, with other codes [561 EZ]- A decrease in 
ion heat flux by a factor of 65% is seen in FigJT3b for the 7^ = 0, R/L Ti = 9.2 case 
between (3 e = — 0.48%. Simultaneously, while the ion heat flux is reduced by f3 e in the 
7# = 0, R/LiTi = 4.6 case, it is not totally quenched. The observation that for /3 e > 
the flux level is diminished over a range of R/Lxi, yet is not totally quenched in the 
vicinity of the ITG threshold for (3 e = 0, is indicative that (3 e > (within the range 
studied) induces a decrease in stiffness as opposed to a threshold shift. Note that the 
results reported in Ref. [55J cannot be compared with those in FigJT3b quantitatively, as 
TEM contributions to the overall turbulence picture may change in particular the (3 e 
dependence of the threshold shift. 

We can thus conclude that electromagnetic effects play a significant role in stiffness 
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reduction for our parameters, even at relatively low values of ft e . While this stiffness 
reduction is not sufficient to fully explain the experimentally observed stiffness reduction, 
it is a factor which must be taken into consideration. 

Additionally, from linear gyrokinetic analysis, (3 e stabilisation of ITG modes has 
been invoked as a possible factor in improved hybrid scenario confinement at ASDEX 
Upgrade and DIII-D, particularly at outer radii (i.e. beyond half radius) [58]. The 
increased non-linear /3 e stabilisation reported here may point to an even greater 
importance of this effect than previously recognised. We note that the /3 e stabilisation is 
expected to be effective up to the recently discovered Non-Zonal Transition (3 e limit |59j , 
beyond which electromagnetic fluctuations effectively short out the zonal flows and lead 
to a significant increase in the saturated level of the ITG turbulent fluxes. This (3 e 
threshold very strongly depends on the background gradients, however, and for typical 
(low) gradients quickly becomes less restrictive than the KBM threshold. Coupled with 
the fact that this effect produces a limit with enormous stiffness, it can therefore be 
expected that standard experimental gradient and (3 e values in outer radii lie below this 
point, putting those cases in the (3 e stabilisation zone. 

In Refs. [53| 155] . an increase of the ratio between the zonal flow shearing rate to 
the unstable mode growth rate (cozf/j) was observed with (3 e . A possible physical 
mechanism for this relative increase in zonal flow activity, based on increased coupling 
between Alfvenic modes and drift waves, has been suggested [60]. In FigJHJwe plot the 
mode amplitude spectra for the 7^ = 0, R/L Ti = 9.2 scan over (3 e shown in FigJTBI 
The amplitude spectra have been normalised to the zonal flow (or rather k y = 0, which 
constitutes a reasonably good measure) amplitudes. Indeed, a relative increase in the 
k y = modes is seen for the electromagnetic cases, which may be related to the ITG 
/3 e stabilisation. Another possible mechanism for increased zonal flow coupling is the 
observed widening of the ITG eigenmode structure observed with increasing (3 e , as shown 
in Fig{T5j The less ballooned structure facilitates the direct coupling to the poloidally 
symmetric zonal modes, similarly to what occurs at low magnetic shear [611 [39] . Further 
work is suggested to shed more light on this topic. 

It is interesting to note that the stabilising effect of flow shear is weakened by finite 
(3 e in the higher R/Lxi case, as seen in FigfT3j In the R/Lxi = 9.2 case, the effect of 
flow shear on the turbulence switches from stabilising to destabilising as (3 e increases. 
However, in the R/L^i = 4.6 case flow shear is always stabilising, and no discernible 
weakening of the stabilisation is seen as f3 e increases. Linearly, the PVG modes are not 
observed to lead to increased stabilisation at increased (3 e . More effort needs to be taken 
in the future to uncover the non-linear effects which either increases PVG destabilisation 
or decreases the ExB stabilisation in the high R/Lti case. 

In Figjl6]we can see, from the entire data set in Fig.l of Ref . [25] . the correlations 
between R/Lxi and (3 e for p = 0.33 and p = 0.64. There is a generally limited but 
positive correlation between (3 e and R/L Ti at p = 0.33, consistent with the reduced 
stiffness. At p = 0.64, (3 e is generally much lower than at p = 0.33. This is expected 
to further increase the stiffness at p — 0.64 as observed experimentally, beyond the 
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Figure 13. Linear (a) and non-linear (b) j3 e scans with 66404 parameters at p = 0.33. 
R/ Lxi and je are varied. Runs were with collisions and numerical geometry. 
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Figure 14. Amplitude spectra from the je = 0, R/Lri = 9.2, non-linear f3 e scan 
displayed in FigfT3b. 



increase solely due to the higher s and q values. We note that both at p = 0.33 and 
p = 0.64, a cluster of five points is visible at the highest respective (3 e values, separate 
from the main trend. These points correspond to discharges with significantly more 
heating power (between 10 — 15 MW of NBI power) and slightly lower magnetic field 
(3 T as opposed to 3.4 T) than the rest of the dataset. The scatter that these points 
induce to the correlation between (3 e and R/Lfi is indicative of the difficulty in making 
a pure comparison of the effect of (3 e throughout the dataset, due to the concomitant 
changes in other plasma parameters and normalized ion heat fluxes. 



Ion temp, profile stiffness: non-lin. gyrokinetic simu. and comp. with exp. 



19 




Figure 15. f3 e scan of ITG eigenmode structure calculated by linear-GENE. R/Lti 
9.2, and ~f E = 0. 
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Figure 16. Correlation between R/Lti and (3 e at p = 0.33 (a) and at p = 0.64 (b) 
from the entire data-set presented in Ref. |25j . 



3. 7. Full simulations including all effects 

Simulations of all four discharges in the data set at p = 0.33 were carried out. We analyse 
the 'high-stiffness-branch' and 'low-stiffness-branch' separately in sections 13.7.11 and 
13.7.21 respectively. Ion heat fluxes from non-linear simulations and experimental power 
balance are compared. The simulations included flow shear, the effect of rotation on 
equilibrium, experimental R/L n , finite /3, collisions, Z e ff > 1, and experimental T e /Tj. 
The effect of Z e ff - which is stabilizing for ITG turbulence - was modelled by lumping 
all impurities into a kinetic fully stripped carbon ion species. The carbon temperature, 
R/Lt, and R/L n were assumed the same as the main deuterium species. Simulations 




Figure 17. Sensitivity of growth rates to Z e f f (a) and T e /Tj (b) from linear Gene runs 
based on 66404 parameters at p = 0.33. Runs were electromagnetic, with collisions, 
and with numerical geometry. 



with varying Z e jf values were carried out, to test the sensitivity of the predictions to 
the uncertainties in the Z e ff profile shape. The growth rate sensitivity to Z e ff and 
T e /Ti for linear Gene runs based on discharge 66404 can be seen in FigJTTl For our 
cases, the Z e ff stabilisation tends to be compensated by the T e /Tj > 1 destabilisation. 
In the non-linear simulations, assuming R/Lt c = for the carbon species instead of 
R/Ltc = R/Lxi altered the bulk ion heat flux by less than 2%. While the impact of 
the fast ions on the Shafranov shift is included in the simulations, fast ion dilution of 
the main ion species was not included. It was judged that the relatively minor (~ 10%) 
impact of this effect on the simulated flux values did not justify the inclusion of the 
fast ions as a separate ion species in the simulations, which significantly slow down the 
calculations. 

3.7.1. Investigation of the low-rotation, high- stiffness branch In section [37TI figure HI it 
is evident that the stiffness of the simulated s/q = 0.7/1.7 curve is less than the apparent 
experimental trend. In this section, we examine the possibility that the higher T e /T, of 
the high flux discharge 73221 in the low rotation branch is responsible for the increased 
flux, through the T e /Ti impact on the ITG critical threshold. It is important to note 
that this significant difference in T e /Tj between the high and low flux discharges in the 
low-rotation branch was not apparent when the data was analysed in Refs. [24, 25J. 
Since that time, the ECE data has been recalibrated. 

A R/LiTcrit oc (1 + Ti/T e ) scaling has been derived both analytically and from linear 
gyrokinetic simulations for the ITG instability [2, 38J. A decreased instability threshold 
leads to increased flux for a given R/Lti value, as long as the stiffness level does not 
change with T e /Tj. It has been predicted by non-linear simulations that the stiffness 
level is not sensitive to T e /Tj within the range relevant for our studied discharges [62J. 
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The simulation results for the 70084 and 73221 discharges are shown in figure [TBI Since 
R/Lti is close to threshold and the transport is relatively stiff, the results are highly 
sensitive to the input parameters. Additionally, the proximity to threshold leads to 
statistical flux variations due to intermittency often higher than the typical 5 — 10% 
level observed for the simulations in this paper. These variations are displayed on the 
plot for these specific cases. For 70084, agreement between the non-linear simulation 
and the experimental observation was found for reasonable departures from the base 
parameters recorded in table El R/Lt% and T e /Tj were both taken at the high end of 
their error bars. Z e ff was taken as 1.4, lower than (Z e ff) = 2.2. This is a reasonable 
assumption since the Z e ff profiles tend to be hollow, and p = 0.33 is relatively close 
to the magnetic axis. Making the same assumptions for 73221 (although maintaining 
the base value of T e /Tj), the simulated flux value was found to be significantly lower 
than the experimental value. Even though T e /Tj is higher in 73221 than in 70084, the 
impact of the higher T e /Tj on the ITG critical threshold is compensated by the lower q 
value calculated by the 73221 CRONOS interpretative simulation compared with 70084. 
However, when increasing the 73221 q value in the simulation to equal the 70084 value 
- an increase of only ~ 15% - the simulated flux value then becomes comparable to the 
experimental value. When assuming the Faraday rotation constrained EFIT g-profile 
for 73211, with s/q = 0.5/1.4, we obtain an intermediate flux level between the 70084 
and 73221 experimental flux values. These tests of the 73221 to the variations of q and 
s constitute a sensitivity analysis of the fluxes to reasonable estimates of the g-profile 
error bar. We thus deem that the T e /Tj increase of the high flux cases in this branch 
compared with the low flux cases is a likely explanation for the seeming anomalously high 
stiffness of this data-set. However, the high sensitivity of the simulated flux - through 
the impact on the critical threshold - to T e /Tj and the g-profile variations within the 
estimated experimental error bars precludes a firm conclusion on this point. The result 
lies within the uncertainties - particularly of the g-profile calculations. In table |3] we 
show the results for all simulations carried out for 70084 and 73221 - beyond those shown 
in Fig. [THl The sensitivities of the flux to input parameters such as s, q, Z e ff, je, and 
R/L n are shown. We note that the results are not highly sensitive to wide variations in 
the R/L n values. 

3.7.2. Investigation of the high-rotation, low- stiffness branch The comparison between 
the Gene non- linear simulations and the experimental heat fluxes for the 'low- stiffness 
branch' is shown in FigJT9j For the high rotation discharges, three separate sets of 
simulations are shown: with the g-profile from the CRONOS interpretative runs and 
Z e ff = 1.9, with the optimistic s/q = 0.2/1.3 assumption and Z e ff = 1.9, and finally 
with the optimistic s/q = 0.2/1.3 assumption and Z e ff = 2.4. The input parameters 
and flux values for these simulations, as well as additional simulations carried out for 
further sensitivity studies and for clarity not shown in FigJT9l are listed in table HI 

For the R/Lti = 6 discharge 66130, the simulation with the base parameters (i.e. 
with the CRONOS s and q values) leads to a flux value x ~ 2.5 above the experimental 
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Figure 18. Comparison of experimental and simulated ion heat flux for discharges 
70084 and 73221 situated on the 'high-stiffness-branch' at p — 0.33 from the dataset of 
Ref.|25j. The 73221 simulation results shown were carried out for three separate s/q 
values to test the sensitivity to the g-profile uncertainties. 



Table 3. Input data and ion heat flux results for discharge 70084 and 73221 non-linear 
simulations. The cases in bold font are the simulations displayed in Fig. 1181 



Shot number 


Zeff 


R/Lm 


R/L„ 


T e /T, 


7E 


s 


g 


qi [gyroBohm units] 


70084 


1.4 


3.5 


1.4 


1.12 


0.07 


0.7 


1.7 





70084 


1.4 


4 


1.4 


1.12 


0.1 


0.2 


1.3 





70084 


1.4 


4 


1.4 


1.12 


0.1 


0.7 


1.7 


8.2 ± 1.4 


70084 


1.4 


4 


1.4 


1.12 


0.07 


0.7 


1.7 


14 ±4 


70084 


1.4 


4 


1.4 


1.08 


0.1 


0.7 


1.7 





70084 


1.9 


4 


1.4 


1.12 


0.07 


0.7 


1.7 





70084 


1.9 


4 


1.4 


1.12 


0.04 


0.7 


1.7 


7.5 ±1.5 


73221 


1.4 


4.2 


2.8 


1.35 


0.02 


0.7 


1.4 


12 ±3 


73221 


1.4 


4.2 


1.0 


1.35 


0.02 


0.7 


1.7 


48 ±2 


73221 


1.4 


4.2 


2.8 


1.35 


0.02 


0.7 


1.7 


40 ±7 


73221 


1.4 


4.2 


3.8 


1.35 


0.02 


0.7 


1.7 


31 ±7 


73221 


1.4 


4.2 


2.8 


1.35 


0.02 


0.5 


1.4 


20 ±2 


73221 


1.9 


4.2 


2.8 


1.35 


0.02 


0.7 


1.4 


1.7 ±0.3 


73221 


1.9 


4.2 


2.8 


1.35 


0.02 


0.7 


1.7 


13 ±3 



Table 4. Input data and ion heat flux results for discharge 66130 and 66404 non-linear 
simulations. The cases in bold font are the simulations displayed in Fig. [TjJ] 



Shot number 


Z 'ff 


R/lm 


T e /Ti 


s 


<J 


qi [gyroBohm units] 


66130 


1.4 


6 


1.25 


0.2 


1.3 


19.4 


66130 


1.4 


6 


1.12 


0.2 


1.3 


12.3 


66130 


1.9 


6 


1.25 


0.7 


1.8 


31.3 


66130 


1.9 


6 


1.25 


0.2 


1.3 


11.1 


66130 


2.4 


6 


1.25 


0.2 


1.3 


7.5 


66404 


1.4 


8.6 


1.14 


0.2 


1.3 


53.2 


66404 


1.9 


8.6 


1.14 


0.4 


1.8 


77.1 


66404 


1.9 


8.6 


1.14 


0.2 


1.3 


33 


66404 


2.4 


8.6 


1.14 


0.4 


1.8 


47 


66404 


2.4 


8.6 


1.14 


0.2 


1.3 


23.8 


66404 


2.4 


7.7 


1.08 


0.2 


1.3 


13.7 
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Figure 19. Comparison of non-linear simulations and experimental results for the 
three separate discharges at p = 0.33. For the high rotation discharges, various sets of 
simulations with varying s, q, and Z e ff assumptions are shown. The Gene simulations 
corresponding to the low and high rotation discharges were carried out with = 0.1 
and 0.3 respectively. 



level. This discrepancy can be reduced by a reasonable variation of input parameters 
around the experimental uncertainties, either for q and s, Z e ff, or R/ Lt%- However, the 
discrepancy between the simulation and the experimental flux for the higher Rj L?i = 8.6 
discharge - 66404 - is significantly greater. For the simulation with the base input 
parameters, the simulated flux is x ~ 5 higher than the experimental value. The 
simulated and experimental flux can only be reconciled by making a highly optimistic 
assumption with regard to the simultaneous variation of R/Lxi, Z e ff, s, q, and T e /Tt 
around their estimated error bars - as seen in the last line of table HJ This is a highly 
unlikely scenario, and leads us to conclude that the stiffness in these specific gyrokinetic 
simulations is higher than the experimental observations for the high-rotation branch. 

The predicted and experimental fluxes for 66404 can however be reconciled by 
both artificially increasing 7^ beyond the measured value from the toroidal rotation, 
and simultaneously ignoring PVG destabilisation. This is shown in an additional set 
of simulations displayed in Figj20l This assumption is consistent with assuming non- 
negligible poloidal rotation. Our original assumption of negligible poloidal rotation 
due to neoclassical damping was justified by nclass [63] neoclassical poloidal rotation 
predictions for the deuterium species within the CRONOS modelling. This is seen in 
Figj21j where the 7^ profiles derived from the nclass predicted poloidal rotation are 
shown. While there is an increase in correlated with increasing R/Lxi as expected, 
the absolute values are - while not entirely negligible for the 66130 and 66404 cases 
- still approximately an order of magnitude below the values necessary to provide 
significant turbulence suppression as observed. However, poloidal rotation values 
significantly above neoclassical values have been observed within internal transport 
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Figure 20. Comparison of flux values from non-linear simulations and experimental 
power balance for the high rotation discharges 66130 and 66404. Here we assumed for 
simplicity T e /Ti — Z e ff — 1. Sets of simulations both including and excluding the 
PVG drive are shown. For each set, additional 66404 simulations with increased 
from 0.3 to 0.6 were carried out. 



barriers (ITBs) [10]. nclass predictions have also been shown to deviate from 
experimentally measured carbon poloidal rotation values at DIII-D [64]. Also at DIII-D, 
the directly measured differences between core main ion and impurity toroidal rotation 
are observed to be inconsistent with the neoclassical predictions in cases with steep 
pressure gradients [65J. For our study, this questions both the validity of the negligible 
poloidal rotation assumption and also the inference of the main ion rotation level from 
the impurity rotation measurements. It is thus of interest to directly measure poloidal 
rotation in the low-stiffness-regime discharges, to examine whether nonetheless any 
anomalous poloidal rotation is observed. It is also of interest to examine a theoretical 
mechanism for generation of anomalous poloidal flow - potentially via a turbulent 
Reynolds stress- particularly for cases with a high degree of external toroidal momentum 
injection. 

The dramatic decrease in flux seen in the 66404 case when transitioning from 
7s = 0.3 to 7s = 0.6 in the no-PVG case in FigEOcan be understood as a consequence 
of the suppression of lower k y modes. In Figj22l a 2D linear R/Lti and k y scan of 
growth rates for discharge 70084 (including realistic parameters) is shown. 70084 is 
shown instead of 66404 for this illustration since the ITG behaviour is similar but the 
TEM driven modes at low R/L^i in the 66404 case complicates the picture without 
changing the qualitative message. No rotation was included in the scan. The contours 
for 7 = (linear threshold), 7 = 0.15 and 7 = 0.3 are highlighted in white. Modes 
with growth rates up to 7 = 0.15, 0.3 can be expected to be suppressed by 7^ ~ 0.3, 0.6 
respectively. Low k y modes drive a significant amount of flux, as represented by the 
7/fc 2 mixing length argument. In Figj22| the estimated low k y cutoff due to a flow shear 
of 7s = 0.6 is thus k y ~ 0.3. In that regime, the flux should be significantly reduced 
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Figure 21. derived from the NCLASS predicted poloidal rotation for deuterium. 
The solid lines are the average values over the 1 s time window studied for each case. 
The dashed lines corresponded to the standard deviation of the profiles around the 
mean during the time window. 

since the high-flux-inducing low k y modes are sheared out of the system. Indeed, when 
comparing the flux spectra maxima for the 'Je = 0.3 and 7# = 0.6 no-PVG non-linear 
cases, the maxima are at k y = 0.25 and k y = 0.35 respectively, in line with the intuition 
gained from the linear runs. 



0.1 0.2 0.3 0.4 




Figure 22. 2D plot of growth rates as a function of R/Lxi (x-axis) and k y (y-axis) for 
discharge 70084 with realistic parameters - identical to the non-linear case in FigfH)l 
The growth rates are in units of [c s /R]. The three white contours correspond to 
7 = 0, 0.15, 0.3 from left to right respectively. 

We now summarise the entire discussion on the low-stiffness question. The predicted 
impact of the differences in parameters between the low and high rotation discharges at 
p = 0.33 were examined in detail with linear and non-linear gyrokinetic simulations to 
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investigate the potential factors leading to the observed reduced stiffness in the high- 
rotation cases. It was found that the differences in R/L n and the effect of rotation 
on the equilibrium have negligible impact on the stiffness for our parameters. The 
effect of rotation itself, and of the fast particle content in the high rotation cases, have 
non-negligible but insufficient impact to explain the observed difference in rotation. The 
impact of q and s on the stiffness level is however significant. The non-linear stabilisation 
of ITG turbulence due to electromagnetic effects (/3 e ) was significant, but by itself still 
insufficient to explain the full stiffness reduction. When self-consistently including all 
effects, the ion heat flux values predicted by the gyrokinetic simulations agreed with the 
observed values in the low rotation case (70084), and were x ~ 2.5 and x ~ 5 higher 
for the high rotation cases 66130 and 66404 respectively. For reasonable variations 
of the input parameters around their uncertainties, the simulated and experimental 
flux values for 66130 could be reconciled. However, for 66404 (which is the highest 
R/L Ti case), the discrepancy between the simulation and the observation could not be 
reconciled by a reasonable variation of the input parameters. Improved agreement for 
66404 could however be obtained by assuming both a 7^ value higher than measured 
from the toroidal rotation profile, and simultaneously suppressing parallel velocity 
gradient destabilisation. These assumptions are consistent with non-negligible poloidal 
rotation, which is an unmeasured quantity for these discharges. Poloidal rotation 
gradients approximately an order of magnitude higher than the predicted neoclassical 
values would be necessary to achieve sufficient impact. Carrying out poloidal rotation 
measurements for this class of discharge is thus recommended, as is the analytical and 
numerical investigation of potentially significant turbulent generation of poloidal flow 
in the relevant parameter regime of these discharges. 

4. Simulations at p = 0.64 

In the previous section, the possible factors leading to a difference in stiffness between 
the low and high rotation discharges at p = 0.33 was investigated. In this section we 
investigate the experimental observation of a lack of stiffness reduction with rotation 
between the classes of discharges at p = 0.64, which attained similar R/Lxi values, 
as seen in Fig{Tb. Non-linear simulations with gene of three of the discharges were 
performed, with parameters matching those at p = 0.64. First, reduced simulations are 
carried out based on 70084 parameters, varying the rotation alone and examining its 
impact on R/Lxi and the stiffness. Then, full simulations are carried out - analogous 
to those in section |3~7T1 - and the Gene predicted ion heat fluxes are compared with the 
experimental values. 

In Figj23] a non-linear R/Lt% scan with various levels of 7^ is shown. The scan 
is based on discharge 70084 parameters, but uses circular geometry, s/q = 2/3, and is 
collisionless and electrostatic. The simulated stiffness is indeed greater than the p = 0.33 
case shown in Figj5]as can be seen in a direct comparison shown in FigfJUfor the 7^ = 
case. Moreover, the degree of experimental 7^ variation between the discharges (between 
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7# = 0.1 — 0.3) is also not sufficient to lead to a significant difference in R/L-n for the 
same level of flux. 

Examining the differences in experimental parameters for all 3 discharges between 
p = 0.33 and p = 0.64 in table [2J we can see that both s and q are higher at p = 0.64, 
and /3 e is lower. All of these differences are expected to lead to higher stiffness in the 
p = 0.64 cases compared with p = 0.33. These qualitative differences in g-profile and 
(3 e between low and high radii are generic (apart from special cases such as in ITB 
discharges), and should hold in general in tokamak discharges. 

In Figj25l the full comparison between the simulations and the experiments is 
shown. These gyrokinetic simulations are electromagnetic, collisional, with numerical 
geometry, include a carbon species at a density consistent with Z e ff = 1.9 for 66130, and 
Z e ff = 2.4 for 70084 and 66404. The simulations include the experimental T e /Tj. For 
all cases, the simulated and experimental ion heat flux agree approximately within 50%. 
This magnitude of difference can be easily reconciled within the reasonable uncertainties 
of input modelling parameters such as R/Ltz, T e /Ti or Z e //, particularly for these stiff 
transport cases. An R/L^e sensitivity check for discharge 66130 was carried out, which 
had the largest relative R/LT e error throughout the data set, as seen in table El It was 
found from the dedicated non-linear simulations that within the possible R/L^e range 
the impact on ion transport is minimal. 



250 




2 4 6 8 10 

R/L 

Tl 

Figure 23. Non-linear R/Lti scan for various levels of je, based on the 70084 
parameters at p = 0.64. Circular geometry, s/q = 2/3, collisionless and electrostatic. 

To summarise, the effect of rotation alone at p = 0.64 is not expected to 
lead to experimentally discernible differences in R/Lt% and stiffness for the range of 
experimental examined. This is in agreement with the experimental trend seen 
in FiglTb. When comparing the full non-linear gyrokinetic ion heat flux predictions 
with the experimental values at p = 0.64, general agreement within reasonable input 
parameter uncertainties is seen for all the discharges, both at high and low rotation. 
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Figure 24. Non-linear R/Lxi scan comparing the stiffness level at p = 0.33 and 
p = 0.64, at 7# = 0, based on the 70084 parameters. Circular geometry, collisionless, 
and electrostatic. 
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Figure 25. Comparison between gyrokinetic simulations and experiment at p = 0.64 
for all three discharges. The experimental values (with the error bars) are shown for 
70084 (red marker), 66130 (green marker) and 66404 (blue marker). The simulated 
values are shown with the same colour coding and marker style for all three discharges. 
Runs were electromagnetic, with collisions, and with numerical geometry. 



5. Conclusions 

Observations at JET have shown evidence of reduced ion temperature profile stiffness 
correlated with low magnetic shear and increased flow shear. The same data-set has also 
raised questions regarding the experimental validation of the Dimits shift paradigm, and 
the low-rotation subset of discharges within this data-set seemed to display higher profile 
stiffness than expected from gyrokinetic simulations. These observations have motivated 
extensive non-linear gyrokinetic simulations to investigate these questions. Simulations 
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using the Gene code were carried out, with parameters based on a subset of these 
JET discharges. Transport sensitivity scans of various parameters that differed between 
the discharges - aside from rotation - were carried out, to assess potential mechanisms 
that may explain the observations. Full simulations including electromagnetic effects, 
numerical geometry, Z e ff, experimental T e /Tj, and rotation were also performed at 
p = 0.33 (in the low stiffness zone) and p = 0.64 for the discharges studied. The 
predictions were compared with the experimental results. The conclusions can be 
summarised as follows: 

(1) For the low rotation cases at p = 0.33, agreement between the gyrokinetic 
simulations and the experimental ion heat fluxes could be obtained within reasonable 
variations of the input parameters within their uncertainties. For the high rotation 
discharges, such agreement could also be obtained for the intermediate Rj ' Lt% = 6 
case, but not for the high R/L Ti = 8 case. Thus, according to the gyrokinetic 
simulations, the measured values of toroidal flow shear are insufficient to explain 
the full transition to the low stiffness regime. While the competition between 
parallel velocity gradient (PVG) destabilisation and ExB stabilisation can reduce 
the stiffness in the vicinity of the turbulence threshold, the predicted flux levels 
themselves are still significantly higher than the experimental values. Improved 
agreement between the simulation and observation for the R/Lxi = 8 case is however 
obtained when both ignoring PVG destabilisation, and assuming ^ E values above the 
measured values. These rotation settings are consistent with non-negligible poloidal 
rotation, which we have assumed to be negligible due to neoclassical damping. It is 
thus of interest to examine both experimentally and theoretically the possibility of 
anomalous poloidal rotation in this regime, potentially due to a turbulent Reynolds 
stress in the presence of external toroidal momentum input. 

(2) The transport sensitivity to R/L n variations, dilution due to fast particles, increased 
Shafranov shift due to suprathermal pressure, and the effect of rotation on the 
equilibrium, were all examined. It was established that none of the above factors are 
sole mechanisms for the transition to the reduced stiffness regime. Their cumulative 
effect is however not negligible - particularly that of fast particles both through 
dilution and an increased Shafranov shift. 

(3) The sensitivity of the transport to f3 e was examined. It was established that 
even for the relatively low j3 e values present in these discharges, the non-linear 
electromagnetic ITG stabilisation is significant. This stabilisation, at least for 
(3 e < 0.48%, is a stiffness reduction as opposed to a threshold shift. The non- 
linear stabilisation is significantly greater than the linear fl e stabilisation, and may 
be related to an increased relative amplitude of zonal modes. Further investigations 
of the parameterisation of this effect is important for incorporation into the 'mixing 
length rule' of quasilinear transport formulations. It is expected that this effect 
would, using such formulations, lead to more optimistic predictions for the energy 
confinement in future devices such as ITER and DEMO, which are not expected to 
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have significant rotation but could still benefit from (3 e stabilisation. 

(4) No clear disagreement is observed between the experimentally observed turbulence 
R/Lti threshold and the upshifted (Dimits shift) non-linear threshold predicted by 
the gyrokinetic simulations. Previously reported results of such a disagreement in 
Refs. [2H [25] were found to be highly sensitive to the precise choice of q values used 
for the simulations. Recently improved data processing methodology has led to a 
revised q value now seemingly pointing to the opposite result of good agreement 
between the experimentally observed and simulated threshold values. However, a 
firm conclusion on this point is not justified considering the sensitivity of the results 
to both q and s. 

(5) For the low-rotation branch at p = 0.33 within the data-set studied, the observation 
of seemingly anomalous high stiffness compared with the gyrokinetic simulations 
is likely explained by a downshift in the ITG critical gradient due to higher T e /T, 
in the high flux cases. However, a firm conclusion in this regard is precluded by 
the high sensitivity of the critical gradient to q and s, and thus to the g-profile 
uncertainties. 

(6) The gyrokinetic predictions and experimental fluxes were also compared at p = 0.64 
for the three discharges. The experimental variation in flow shear between the 
discharges was not predicted to be sufficient to lead to a discernible difference in 
R/Lxi - in agreement with the observations. The simulated and experimental ion 
heat fluxes for all examined discharges all agreed to within approximately 50%. 
This degree of discrepancy can be explained by reasonable variations of the input 
parameters within the experimental uncertainties. 
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